Exercise: ETA vs Manual assembly

Aim of this exercise is to solve a Laplacian problem first with ETA and then by writing a custom assembler class, which does not use ETA for assembling the stiffness matrix. This is useful because, while ETA is a powerful tool for expressing directly the weak formulation of the problem we want to solve, minimizing the time spent from the mathematical formulation to the its implementation, there are situations for which maximum performance can be achieved only by directly coding the assembly loops.

The problem to solve is a Poisson problem with a specified RHS and homogeneous essential boundary conditions. The source term at RHS and the analytical solution can be found below. In order to assess the correctness of the code, we will check the error in L^2 norm between the obtained numerical solution and the analytical solution.

Before doing this exercise, please follow all the parts of the tutorial, as we heavily rely on them for the whole exercise.

This is the list of things to do:

Setup the project as an external application using an installation of LifeV. To this aim, into a new folder create a CMakeLists.txt file which:
- finds LifeV installed in a certain directory;
- compiles the files main.cpp and MyAssembler.cpp and creates a test executable.
See this as a reference.
By following the relevant tutorial parts, start coding the main source file main.cpp:
- add the base necessary LifeV headers
- setup MPI and the communicator
- load a datafile from an input file (e.g. input.in)
- load and partition the mesh cube.mesh
- create the Finite Element spaces
- allocate the system matrix
Implement the assembly of the stiffness matrix based on ETA.
Assemble the RHS.
Apply the boundary conditions to the boundary with flag 3.
Solve the linear system.
Compute the error in L^2 norm against the analytical solution (see below)
Export the solution and the analytical solution.
Visualize the solution with Paraview.

Your solver should now be up and running.

Second part of the exercise:

Add a header file MyAssembler.hpp hosting the definition of the custom assembler class. See the template file provided.
Add the implementation of the custom assembler class in a file named MyAssembler.cpp. See the template file provided.
Continue the main.cpp with the code performing the assembly:
- instantiate a MyAssembler class
- call the prepare function
- call the assemble function
Compare the solution and the errors obtained with the two assembly methods.

Tips and references

In order to do the exercise, reading the LifeV documentation can be beneficial. In particular, we highly suggest to read the following references:

CurrentFE: class holding the element currently being processed, providing methods to get the value and derivatives of basis functions in the quadrature points, amount the others (reference). In particular, read the role of the flags in the update process, as detailed in the section Detailed description.
MatrixElemental: class holding the element matrix (reference).
Reference for pushing local element matrices to the global stiffness matrix: Assembly.hpp.

Moreover, please note the following suggestions:

To update the current element when the e-th element is processed: call currentFE->update(...) with the element objtained by calling fespace->mesh()->element( e ) and the update flags set accordingly (see the documentation)
To reset a matrix: matrixPtr->zero()
To get the k-th derivative of the i-th basis function in the q-th quadrature point: currentFE->dphi( i, k, q )
To get the (weighted) determinant of the Jacobian in the q-th quadrature point: currentFE->wDetJacobian( q )
To set the value val into the entry (m, n) of a matrix: matrix(m, n) = val.
This particular problem is such that each component of the unknown is independent from the others. Therefore, it is possible to assemble the local element matrix for just one component and then copy it for the other components. In order to do this, we suggest to assemble the local element matrix for one component into the provided M_tmpMatrix and then copy it into M_localMatrix into the right blocks. This is the code to copy the temporary matrix into the block (m, n) of the local matrix:
```
#!C++
    // 0 <= m, n <= fieldDim - 1
    MatrixElemental::matrix_view mat = M_localMatrix->block( m, n );
    mat += (*M_tmpMatrix);
```

To push the current local element matrix to the global stiffness matrix: see the documentation of the assembleMatrix(...) function defined in Assembly.hpp. In particular, for this problem you can use the following code:

#!C++
    // f = index of the component of the field being considered (0 <= f <= fieldDim - 1)
    assembleMatrix( *matrix, *M_localMatrix, *M_currentFE, *M_currentFE,
                    M_fespace->dof(), M_fespace->dof(),
                    f, f, f * numTotalDofs, f * numTotalDofs );

Computing the error

Steps to do in order to compute the error between the obtained numerical solution and the analytical solution:

Interpolate the exact solution on the FE space (see the tutorial on boundary conditions).

Compute the (local) error by performing a numerical integration on the (local) domain using ETA:

Real errorL2 = 0.0;
{
    using namespace ExpressionAssembly;
    integrate ( elements ( local_mesh ),
                uFESpace->qr(),
                ...
              ) >> errorL2;
}

Accumulate the local errors from all processes to compute the global error:

Real totalErrorL2 = 0.0;
// let all processes arrive to this point
Comm->Barrier();
// gather the local errors and sum them up
Comm->SumAll (&errorL2, &totalErrorL2, 1);

RHS and analytical solution

// To define the source term and the exact solution, we exploit their separability
// Function defined in the xy plane (w can be x or y)
Real funXY(Real w)
{
    const Real L2( 1.5 * 1.5 );
    return 0.5 * (L2 - w * w);
}

// function defined along the z axis
Real funZ(Real z)
{
    const Real L( 3.0 );
    return 0.5 * ( L - z ) * z;
}

// Raw function that describes the source term
Real fSourceTerm ( Real /* t */, Real x, Real y, Real z, ID /* i */ )
{
    return funXY(y) * funZ(z) + funXY(x) * funZ(z) + funXY(x) * funXY(y);
}

// Exact solution of the Laplace problem
Real exactSolution( Real /*t*/, Real x, Real y, Real z, ID /*i*/ )
{
    return funXY(x) * funXY(y) * funZ(z);
}

MyAssembler.hpp template

#!C++

#ifndef H_MYASSEMBLER
#define H_MYASSEMBLER

#include <lifev/core/mesh/RegionMesh.hpp>
#include <lifev/core/fem/FESpace.hpp>
#include <lifev/core/fem/Assembly.hpp>

namespace LifeV
{

class MyAssembler
{
public:
    typedef RegionMesh<LinearTetra> mesh_Type;
    typedef MapEpetra map_Type;
    typedef MatrixEpetra<Real> matrix_Type;
    typedef std::shared_ptr<matrix_Type> matrixPtr_Type;
    typedef FESpace<mesh_Type, map_Type> fespace_Type;
    typedef std::shared_ptr<fespace_Type> fespacePtr_Type;

    void prepare( fespacePtr_Type fespace );

    void assemble( const matrixPtr_Type &matrix );

protected:
    fespacePtr_Type M_fespace;

    std::unique_ptr<CurrentFE> M_currentFE;

    std::unique_ptr<MatrixElemental> M_localMatrix;

    std::unique_ptr<MatrixElemental::matrix_type> M_tmpMatrix;
};

} // namespace LifeV

#endif // H_MYASSEMBLER

MyAssembler.cpp template

#!C++

#include "MyAssembler.hpp"

using namespace LifeV;

void MyAssembler::prepare( fespacePtr_Type fespace )
{
    M_fespace = fespace;

    M_currentFE.reset( new CurrentFE( M_fespace->refFE(), M_fespace->fe().geoMap(), M_fespace->qr() ) );

    UInt nbFEDof = M_currentFE->nbFEDof();

    M_localMatrix.reset( new MatrixElemental( nbFEDof, M_fespace->fieldDim(), M_fespace->fieldDim() ) );    

    M_tmpMatrix.reset( new MatrixElemental::matrix_type( nbFEDof, nbFEDof ) );
}

void MyAssembler::assemble( const matrixPtr_Type &matrix )
{
    const UInt numElements = M_fespace->mesh()->numElements();
    const UInt fieldDim = M_fespace->fieldDim();
    const UInt numTotalDofs = M_fespace->dof().numTotalDof();
    const UInt numLocalCoor = M_currentFE->nbLocalCoor();
    const UInt numDofsPerElem = M_currentFE->nbFEDof();
    const UInt numQuadPt = M_currentFE->nbQuadPt();

    // TODO: write your code here

}

Wiki

lifev-release / tutorial / Manual assembly